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^ ' The dynamics of pliase field crystal (PFC) modeling is derived from dynamical density functional 

theory (DDFT), for both single-component and binary systems. The derivation is based on a trunca- 
tion up to the three-point direct correlation functions in DDFT, and the lowest order approximation 
using scale analysis. The complete amplitude equation formalism for binary PFC is developed to 
describe the coupled dynamics of slowly varying complex amplitudes of structural profile, zeroth- 
O ■ mode average atomic density, and system concentration field. Effects of noise (corresponding to 

' stochastic amplitude equations) and species-dependent atomic mobilities are also incorporated in 

this formalism. Results of a sample application to the study of surface segregation and interface 
intermixing in alloy heterostructures and strained layer growth are presented, showing the effects of 
different atomic sizes and mobilities of alloy components. A phenomenon of composition overshoot- 
ing at the interface is found, which can be connected to the surface segregation and enrichment of 
, one of the atomic components observed in recent experiments of alloying heterostructures. 
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I. INTRODUCTION 



Understanding the formation and dynamics of complex spatial structures or patterns has been of continuing interest 
^ . due to the fundamental importance of predicting and controlling system properties and material functions. However, a 
^ ■ comprehensive understanding is hindered by the fact that the processes involved are usually nonlinear, nonequilibrium, 
. can span a variety of length and time scales, and are highly influenced by the complex coupling with materials 
' growth and processing conditions. Typical examples include the growth of strained solid films and the formation 
of nanostructures such as quantum dots or nanowires, which involve the interplay among microscopic crystalline 
• ■ structure, mesoscopic or nanoscale surface pattern, topological defects (e.g., dislocations), as well as various growth 
[ parameters such as temperature, misfit strain, growth rate, and film thickness [H-Q- The system dynamics and 
evolution are further complicated in alloy samples, due to the additional coupling to spatial/temporal variation of 
alloy composition particularly in the case of phase separation ^3, ^ . 

To address these complex phenomena a variety of theoretical modeling and simulation methods have been developed, 
which can be roughly characterized via the level of description that they focus on. At the microscopic level capturing 
crystalline details, atomistic modeling techniques such as Monte Carlo (MC) or molecular dynamics (MD) have been 
5J] widely adopted. For example, nanostructure (e.g., islands/pits) formation during strained film epitaxy has been 
5^ ] studied via the kinetic MC method incorporating elastic interaction and strain energy , while detailed structure 
and dynamics of crystal defects like grain boundaries and dislocations have been simulated by MD [1, [l^]. However, 
the limitation of small length and time scales addressed in these atomistic methods leads to large computational 
demands and hence the restriction of system size and evolution time range that can be accessed. Such limitation can 
be overcome via continuum modeling methods, including continuum elasticity theory used in strained film growth 
[sl. [Tll - [l7j and the well-known phase field models that have been applied to a wide range of areas such as crystal growth, 
nucleation, phase separation, solidification, defect dynamics, etc. |18l - [23 |. These continuum approaches feature coarse- 
grained, long-wavelength scales and diffusive time dynamics, but are not formulated for the short-wavelength scales 
associated with microscopic crystalline details. 

To incorporate the advantages of these approaches and hence be able to simultaneously model crystalline details on 
length and time scales of experimental relevance, the phase field crystal (PFC) 23 25] model and the related amplitude 
representation [2^ - l30j were developed recently. The PFC model incorporates the small length scales of the crystalline 
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structure with diffusive time scales by describing tfie dynamics of tfie atomic number density field p, a continuum field 
variable that is spatially periodic at atomic length scales in crystalline state [H,[23| • To alleviate the limitation imposed 
by the necessity of describing small length scales, an amplitude representation was developed to describe slowly varying 
envelope or amplitude functions while maintaining the basic features of crystalline states, particularly elasticity, 
plasticity and multiple crystal orientations. Both the original PFC and corresponding amplitude representation have 
be extended to binary alloys 



phase of the density field 
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l32| . In the binary case the amplitude representation describes the amplitude and 
30| and also the concentration profile (3^ , which is assumed to vary on "slow" scales 
compared to atomic lattice spacing. A wide range of phenomena has been studied via this PFC method for both pure 
and binary material systems, including solidification [25l |3 2|. [331 . g rain nucleation and growth 26|, 30l . [s^ . l35|. phase 



segregation [25|, , quantum dot growth during epitaxy [3|, |32, |36|, |37| , surface energy anisotropy |38l, 133] , formation 
and mel ting of dislocations and grain boundaries |4CH43| , commensurate/incommensurate transitions [44, 45], sliding 
friction [46l. l47j . and glass formation [48]. In addition, recent work has been conducted to extend the modeling to 
incorporate faster time scales associated with mechanical relaxation [igllsoj . and to develop new efficient computational 
methods [2i,|5lH5|. 

The PFC model can be connected to microscopic theory via classical density functional theory (DFT) of freezing 
[25l . [ssl . [53 - [56| . It has been found that the PFC free energy functional can be derived from classical DFT for either 
pure materials or binary mixtures [57l - l62| , by approximating the two-point direct correlation function with a truncated 
Fourier series and expanding the ideal-gas part of the DFT free energy functional in a power series of p and ip (up to 
4th order) 25]. While this connection provides insight into the parameters that enter PFC models, the approximations 
used are quite drastic and the resulting model is a poor approximation of classical DFT [56| . A similar connection 
could be made with the atomic density formulation of Jin and Khachaturyan jpsj which is similar in form to the 
classical DFT, although the parameters that enter are given a different physical interpretation. 

The main difficulty in directly simulating classical DFT is that the solutions for p are very sharply peaked around 
the lattice positions (at least in metallic crystals), while simple PFC models predict very smooth sinusoidal profiles. 
This difference makes numerical simulations of a simple PFC model much simpler than classical DFT as the former 
model's grid spacing can be a factor of ten larger than the latter's, so that in three dimensions a PFC model can 
simulate systems three orders of magnitudes larger than classical DFT with the same memory requirements. In 
addition it has been shown that a simple PFC can be adjusted to match many material properties, such as surface 
energy and its anisotropy, bulk moduli, and the miscibility gap in three-dimensional (3D) bcc iron [56j and the velocity 
of liquid/solid fronts in two-dimensional (2D) hexagonal crystal of colloids 33]. 

Another benefit of PFC modeling is the ability to efficiently simulate microstructure dynamics. At present, PFC 
dynamics has been largely introduced phenomenologically using time-dependent Ginzburg-Landau type dynamics 
[2^ [25j . Recent progress includes the derivation of hydrodynamic evolution equations for crystalline solids based 
on the Poisson bracket formalism and the simplification to PFC equations [HJ, [5^ . Very recently research has been 
conducted to connect the PFC-type models with microscopic dynamics (Smoluchowski equation) via dynamical density 
functional theory (DDFT) These results were also based on the truncation of DFT free energy up to two-point 
correlation function, and for single-component systems. 

In this paper we provide a systematic derivation of PFC dynamics from DDFT, for both single-component and 
binary systems that involve the evolution of atomic number density and alloy concentration fields (see Sec. |lT|. Our 
derivation includes contributions from three-point direct correlation functions, which have been shown important for 
the DFT calculations [64, 6^]. The original PFC models can be recovered via the lowest order approximation of our 
DDFT results, with the PFC parameters being connected to quantities of DFT correlation functions. Our calculations 
can be directly extended to incorporate fourth and higher order correlation functions in DFT. 

To complete the PFC methodology for binary systems, the full amplitude equation formalism is established for a 
2D system with hexagonal/triangular crystalline symmetry. It incorporates the effects of species-dependent atomic 
mobility and average (zeroth-mode) atomic density that are usually coupled with the dynamics of structural amplitudes 
and concentration field during system evolution but absent in previous studies of binary PFC. As shown in Sec. IIIIl the 
standard multiple-scale expansion is first applied to derive the lowest order amplitude equations, followed by a hybrid 
approach that we develop here to obtain the equations incorporating all orders of expansion. Furthermore, stochastic 
amplitude equations are derived for both single-component and binary PFC models, showing the corresponding noise 
dynamics as well as its coupling due to different atomic mobilities of system components (see Sec. lIVp . 

As has been discussed in previous research, the advantage of the amplitude equation representation can be revealed 
via its large increase of computational efficiency due to the large length and time scales involved [1^, [13, [s^ and also 
its amenability to advanced numerical schemes such as adaptive mesh refinement method [28j . Furthermore, these 
amplitude equations are more amenable to analytic calculations as shown in recent studies of surface nanostructure 
formation in strained epitaxial films [1, [13] as well as in most recent results for establishing the correspondence 
between PFC type models and traditional phase field approaches [12] • To further illustrate these advantages, in Sec. 
IVl we present a sample application of the derived binary PFC amplitude equations to the phenomenon of surface 
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segregation and alloy intermixing. This is of particular importance in material growth (e.g., group VI or III-V 
semiconductor thin film epitaxy |66l473j ). but rather limited information and understanding is available to date. We 
focus on both liquid-solid(crystal) coexistence profile and the coherent growth of strained solid layers, and show the 
control of intra- and inter-layer diffusion by varying material parameters including solute expansion coefficient (due 
to different atomic sizes), misfit strain in alloy layers, and the mobility difference between alloy components. This 
study provides an understanding of mass transport mechanisms during material growth and evolution. The dynamic 
processes of strained layer growth as well as the associated composition overshooting phenomenon are obtained in 
our calculations in Sec. |Vl The results are compared to experimental findings of vertical composition segregation 
or surface enrichment as widely encountered during the growth of various alloy heterostructure systems such as 
InAs/GaAs, Ge/Si, GaAs/GaSb, InP/InGaAs, etc. 



II. DERIVATION OF PFC DYNAMICS VIA DYNAMICAL DENSITY FUNCTIONAL THEORY 



A. Single-component systems 

We start from the DDFT equation governing the evolution of a time-dependent local atomic number density field 



dt 



5J~ 

Mp{r,t)V — 



(1) 



which was first proposed phenomenologically [TJ, UM ^^'^ later derived by various groups via microscopic Brownian 
dynamics [76l - l78| and Hamiltonian dynamics and hydrodynamics [t^ (see Ref. [s^ for a brief review) . The DDFT 
equations for binary A/B systems are similar to Eq. ([!]), with p{r,t) replaced by Pi{r,t) {i — A, B; see Sec. IIIBI 
below), which has also been derived recently from Brownian dynamics (the Smoluchowski equation) [80|. In Eq. ([T]) 
the mobility is M = D/ksT, where D is the diffusion coefficient and T is temperature. 
In classical DFT the free energy functional J" can be expanded as [13, Ull 



j dr[pln{p/pi) - 5p] - j dridr25p{ri)C^'^\ri,r2)5p{r2) 

^ / dridr2dr:iC^^\ri,r2,r:i)5p{ri)5p{r2)5p{r:i) + • • • , (2) 



where pi is the reference liquid state density taken at liquid/solid coexistence, Sp = p— pi and C^") is the n-point direct 
correlation function of the liquid phase at pi . It is important that the correlation functions are taken from the liquid 
state to maintain rotational invariance. Details of the correlation functions depend on the specific material systems 
studied and are usually calculated via various approximations [s^, |6ll, [g^ . Following the original PFC approach [l^ , 
the Fourier component of the two-point correlation function C^^^ is expanded as a power series of wavenumber q to 
fit up to its first peak, i.e., 

(7(2) (q) = -(7o + (72^2 - (74^4 + ... , (3) 

where Cq, C2, and Ci are fitting parameters that can be connected to material properties such as isothermal com- 
pressibility of liquid phase, bulk modulus and lattice constant of crystal state [2^, |5^. For the three-point correlation 
function C3 , its Fourier transform yields 

C(3Hri,r2,r3) = ^^ j dqdq'e'i<^'—'h'i<^'-'^'')C^^\q,q'). 

The simplest approximation is to keep only the zero wavenumber mode, i.e., 

C('Hq,q')^C^'\q = q' = 0)^-Ci'\ (4) 



as adopted in the DFT studies of hard spheres [6J, |8]| and Lennard- Jones mixtures [82j. This can be justified from 
the previous results of hard-spheres DFT calculations that nonzero wavenumber components of C^^^ have been shown 
to yield minor contributions |63 . [8l| and that as order n increases, the oscillation details of (7(") become less and less 
relevant compared to the zero wavenumber mode [ssj l. 
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Defining the rescaled atomic density field n = ip — pi)/pi and using the approximations ([3]) and ([4]), the free energy 
functional 1^ becomes 



AJ^/pikeT = j dr 
where AJ^ = F[p] — F[pi\, and 



1 



(1 + n) ln(l + n) + ^B^n {2R'^V'^ + iJ^'V*^ 



1 



B', 



B'l = piCo ^B'~l, 



B^ 



piCll^Ci, 



R 



B = pfcl,'^/2 



Substituting Eq. ^ into the DDFT equation ([T]), which can be rewritten as 



dn 
'dt 



= M'V 



{l + n)V 



ST 



(with M' = M/pi),we find [8: 
dn 



(B^ - B^)n + B^ {R^V^ + lfn + rn^ + vn^ 



nV{R^W^ + lfn } 



(5) 



(6) 



(7) 



(8) 



where r = -(S^ - B^ + l)/2 + B, v = 2B/3, and we have used the relation M = D/ksT. Note that if only the 
two-point correlation function in the DFT free ener gy (Ht was used it would yield B = v = 0, and Eq. ^ reduces to 
a form equivalent to the PFCl model given in Ref. 33|. However, this would then be a 2nd-order dynamic equation 
due to the absence of term, and as found in our numerical tests, is more difficult to converge at long enough time 
compared to the full 3rd-order Eq. 

It is convenient to rescale Eq. via a length scale R, a time scale R^/DB^, and n — > yJvjB^ n, leading to 



dn^ 
dt 



-en + (V^ +qofn + g2n^ + 



5oV 



where 



e={B- -B')/B-, 90 = 1, 52 = r/V^, go = /S^. 



(9) 



(10) 



The original PFC equation is recovered by considering that (V • [nV(V^ + ^o)^"^]) of higher order compared to 
term V^( V^ This can be obtained via a simple scale analysis: n ~ 0{e^^^) and (V^ + 9o)^" ~ 0{e^^^) (see 

also Sec. IIII Al for more detail of scale expansion). Thus to lowest order approximation, Eq. ^ can be reduced to 
the original PFC model equation that has been widely used: 



dn^ 
dt 



+ (v2 + go)^" + .92"-^ + "-' 



(11) 



This derivation procedure can be readily extended to incorporate higher order direct correlation functions of DFT 
(e.g., four-point, five-point, etc.) and thus to include higher order terms such as n"*, n^, in the PFC model. Similarly, 
these high-order correlation functions can be effectively approximated to lowest order via the zero wavenumber modes, 
based on the recent DFT calculations [s^. For example, the contribution T^'^'> of free energy functional from the 
four-point correlation function is given by 



T^^'^ /ksT = j dridr2dr3driC''^'^ {ri,r2, rg, r4)Sp{ri)Sp{r2)Sp{r3)Sp{r4). 
Assuming C^^^q, q', q") ~ C^'^\q = q' = q" = 0) = —C'^\ the free energy functional ([5]) is now 
AT/pikeT ^ I dr 



(1 + n) ln(l + n) + ]^B^n {2R^V^ + R^V'^) n + ]^B[n^ + ^Bn^ + ^B^n'^ 



where -B4 = pfcjf^ /6. The dynamic equation for n would then be 



^ = Dh^ [-(S" - B')n + B^ (i?2v2 + lYn + Tn^ +vn 
dt I I ^ ' 



^ ' un^ 



{R'^V^ + lfVn I 



(12) 



(13) 



(14) 



where v = 2_B/3 -I- B4 and u ~ 3B4/A. Again the last term of Eq. is of higher order and can be neglected in the 
lowest order approximation. 
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B. Binary systems 



For a binary system with components A and B, the DDFT equations describing the dynamics of the respective 
atomic density fields pA and ps are given by (80| 



dpA 

dt 



MapaV 



ST 

SpA 



dpE 
dt 



= V • 



MbPbV 



SpB 



(15) 



where Ma(b} is the atomic mobihty for specie A (B). The corresponding classical density functional free energy 
(hereafter referred to as "DFT" for short) is of the form [59l - l62| 



Pt In ^ - Spi 
Pi 



F/kBT ^ f dr 

i=A,E 

-Y—. I dri - ■ -drr, ^ C'i".j('^i> ' ' ' > r„)(5p,(ri) • • • (5yOj(r„), 



(16) 



i,...,j=A,B 



where p] is the reference liquid state density of component i, Spi = pi — p], and C^"'^j refers to the n-point direct 
correlation function between components i, j = A, B. Up to three-point correlation functions, we have 



TjkBT = J dr [pA\n{pA/pf) ~ 6pA + PB^n{pB/pf) ~ Spb] 
dridr2 



(2) (2) (2) 

SpA{ri)Cxx{ri,r2)SpA{r2) + ^PB(ri)C^|;(ri, r2)(5ps(r2) + 2SpA{ri)C\^{ri,r2)SpB{r2) 



CAAAi'^i^'^2,r3)SpA{ri)SpAir2)SpA{r3) + C''g!gg{ri,r2,r3)5pB{ri)5pB{r2)5pB{r3. 

(3) 



(3) 



- j dridr2dr3 

+3C!L4B(''i'''2,r3)(5pA(r-i)^PA(r2)(5pB(T-3) + 3C)^^^{ri,r2,r3)6pA{ri)SpBir2)SpBir3) 



(17) 



Similar to the single-component case discussed in Sec. Ill A[ the correlation functions (ri,r2) and Cj-^^(ri, r2, ra) 
(«, J, k — A, B) are expanded in Fourier space as 



(3), 



ijk 



ijk\ 



-a 



ijk 



As in the original binary PFC model, we introduce an atomic density field n and a concentration field if) via 



p- pi PA+ Pb - Pi 



Pi 



Pi 



Pa - Pb pa - Pb 



P 



PA+ PB 



where pi 



pf , and hence 



PA = |(l + n)(l + V^), p5^|(l + n)(l-V'). 



(18) 



(19) 



(20) 



Substituting Eqs. (fT8|) - (l20l) into (|17p . we can express the free energy functional in terms of n and ip: 



AT/pikBT = 



dr I (1 + n) ln(l + n) + -(1 + n) [(1 + ^p) ln(l + ^p) + {I - ijj) ln(l - -0)] 
+ (i±Z^ {2B%m^V^ + B^im'^') " + f |V[(1 + nW + I (V2[(l + n)0])^ 



(21) 



where 
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EL 

4 



5Co 



Pi - Pi 



(3) 



P2 + Y^dGo 



B[{i>) = B^(^) - 1 - i?^ - 1 + B^tA + ^ ^^^3 



Pi 

fH 
4 

d 
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Pi ~Pl 



fH 
2 



.(3) 



Pi - Pi 



~(3) 



AC, 



(3) 



PiSC'o 



(3) 



(3),/,3 



^(3)j,3 



ACn 



Pi - Pi 



AC, 



(3) 



A/3 = /3i-/32 = ^<5Cf , 



B2 = 3A/3, 



2/3o + ^Bi, 



B: 



Pi C2+<5C2^/2 



2 - 4/3i - 3/32, 
2 



4(c4 + (5C4V'/2 



-2 
P;C'2 

4C4 



1 



.C2 



= 2/33, 
2C4 y 



\ 



2(c4 + ^C4V'/2 



is: 



^(C2 



C2 + 5C2^/2 
1 



2C4 
\ C2 



-5C2V^) =i?-i?o^(H-a2V'), 



Bi?i?;5(l + a4V') 



1 / <5C4 ^^2 \ , _^ 

a2 = <5C2/2C2, 
Q!4 = 5Ci/2Ci^ 



— Rn 



Pi 



AC2 



Pi A,/^ 

K — — AC4. 

4 



(22) 



In the above formulae, the following has been defined from the correlation functions: 



C 

C(3) . 



/^(2) , ^(2) , r,^(2) \ 



/^(3) 1 /^V' 



-(3) 



(5C = C 

.(3) 



:2) 



C 



(2) 
SB' 



AC = C 



(2) 
AA 



C 



(2) 
BB 



2C 



r'(3) _ r^i^) 

^ — ^AAA ~ ^BBB 



^(3) 



c 



(3) 
AAB 



C 



(3) 

ABBi 



(23) 



jc(3) = ci^l^ + cL'^ 



'BBB 



■'AAB 



.(3) 
''ABB' 



— ^ AAA R R R •J"-^ A A-R ^ 'JL' 



'BBB 



-(3) 
''AAB 



(3) 

ABB' 



and the " " " in Eq. (1^^ refer to the Fourier coefficients in the expansions of Eq. (|18l) where the numerical subscripts 
on the coefficients refer to the appropriate power of the expansion. For binary alloys the lattice constant is often 
approximated by Vegard's law, i.e., i? ~ i?o + -RiV' — Ro{^ + aip)- In this expansion, near ij; = the solute expansion 
coefficient a is expressed as 



a = Ri/Ro = -{04 - 0:2). 



(24) 



(In the dilute limit (i.e., tp ±1) it would be simple to expand R around ^p ~ ±1 to obtain the solute expansion 
coefficient as well.) Using the simplification adopted in the original binary PFC [s^, it is assumed that ~ Bq and 
R^ ~ ^0(1 + 2(3!^/;), -R'* — ^0(1 + 4:Ctip) via expansion, which corresponds to the assumption of a2 — 2a and ~ 4a 
as obtained from Eqs. and ((211) . 

In terms of the above definitions, the time derivatives of the variables n and ijj defined in Eq. ()19p are given by 



dn 
'dt 



1 f dpA_ ^ dp_ 
dt 



Pi 



dt 



dip _ 1 
'dt ~ piil + n) 



(l-^)^-(l+v^)^ 



(25) 



From the DDFT equations (jlSp . the dynamics for n and thus become 

dn/dt^ M1V1+M2V2, 

d^j/dt = [{M2 - Miij)Vi + {Ml - M2ip)V2] 



(26) 



where 



Ml = ^kBT (Ma + Mb) , M2 = ifc^T {Ma - Mb) , 



(27) 



and 



2?i = 



1 



PikBT 
1 

PikBT 



V • 

V • 



(1 + ")V§^ 

on 

(1 + n)V'V — 
on 



+ V 



(l+n)(l-l/^2)y 



1 

1 + n S^p 



(28) 



Using the free energy functional (|2T|) as well as Eq. (|28t , the results of 2?i and 2?2 (keeping all the terms) are 
-(So" -B'o + l)/2 + Bo + f Bf/2 + Bi) V + + V^'^ 



+ n)V+/?o^ + -(/3i+A/3)V' 



+i?o" v2 + 1)' n + Bo" (a2i?gv2 + ^i?^V4) [(1 + n)i;] } 

+V • [nV [b^ (i?2v2 + 1)2 n + Bo" («2i?^V2 + ^i?4v*) [(1 + n)^ } 

+V • {(1 + n)7AV [Bo" (a2i?gv2 + yi?^V^) n + {-K\/^ + ^tV^) [(1 + n)^]] } , 



(29) 



/?o« + (/3o/2 + B1/3) n2 + ^Bm^ + (1 + /32 + 2A/3n)(l + n)V' 



-B, 



("2^?, 



2v72 , "4 „4v74 



+V • {nV [(1 + n)iP2 + 2Apn)^p + /33(1 + n)^^,^ + Bo" (aziJgV^ + y i?^V' 



-i^V^ + kV^) [(1 + n)tP]] } + V • jv'V 



(1 + n)-iJjV 



2 ~ 

-Binip 
2 ~ 



(BJ - B* + l)n + Bo + -Bi^ + BJ (i?f5V^ + l) n 



BS[a2RlV' + ^Rtw') [{l + n)^} 



(30) 



At this point in the derivation it should be noted that no additional approximations beyond those going into the 
expansions of Eq. ^TE\\ have been introduced. 

1. Non-dimensional form of model 

To simplify the results, the above binary PFC equations can be rescaled via defining a length scale i?o, a time scale 
Rq/MiBq, n y^v/B^ n, and V — > ^J'^l^t V': yielding 



dn/dt = Vi+mV2, d^j/dt^ 



1 



1 + gon 



[{m - goip)'Di + (1 - mgoip)'D2] , 



(31) 



where 



M2 _Ma-Mb _ [be _2~_pf ^(3) 
Ml Ma + Mb \ v 6 3 



(32) 
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If keeping only terms up to 3rd order quantities of n and ■0, the results of Vi and V2 are rescaled as 



I -en + (V^ + qlfn+ {gitjj + gi}?) n + (52 + 92^) + 
+# + v,i^' + u,i;' + (a^V^ + ^V^) [(1 + 5on)V] + 50^ {a^^^ + y V^) n] 
+goV • {nV [(V2 + qlf n + (a2V2 + ^V^) ((1 + 50^)^)] } 
+.goV • [tAV (-/^oV + KoV*) ((1 + gon)4,)\ - goV ■ [{Vyj) (aaV^ + 



(33) 



|.gr7. + (1 + gan) + y V"*) n + {2vi'ip + W21I?) n + {v2 + gtp) r? + 

+woV + t;oV'^ + uoV'^ + (-i^oV^ + KoV) [(1 + .gon)V]} 
+30V • {nV [wnV' + (-J^oV^ + KoV*) ((1 + g^n)^,)] } - .go V 

+50V • {V'V [(-e + 717/') n + 7271^ + (V^ + g^)^ n + [a^^'^ + yV^) ((1 + 50")^-)] } 



(V7^) (a2V2 + ^V* ) n 



where the rescaled parameters are 



(V2 + g2)^n+(a2V2 + ^V4) 7/]} 



(34) 



Wo 



Bl ' 

^0 



go = 1, 3 



.92 



go ( 2^0-1 



30, 



"0 = J<o 



K 



9u 



m n2 



AC2 
2^2 



Kq 



31 



Bi 



92 



Bi + 2Bi 



W2 



+ 2^1 
2i?o^ 



V2 



2v 

'^o/2+_Bi/3^ 



50, 



Ui 



93 



2Bi 



90, 



71 = 3(53 - ui)/go, 



9v ' 

72 = 32 - Vi. 



Note that from Eqs. 



and (dH]), Bq can be rewritten as 



~ 2 

4^4 



pi {Co + C 

max I 7 



(35) 



(36) 



where C,nax is the maximum of the first peak of the two-point correlation function C in Fourier space. If |Ap/ 

- (3) 

Ipi^ - pf I « \Co/Co \,B'o^l + piCo from Eq. 



^0 ~ Bp 



and thus 



(Co 



Co + Cn 



(37) 



Usually Cniax <C Cq. particularly when close to the melting point T,„, and hence e can be viewed as a small variable 
(also used in amplitude equation expansion given below), proportional to (T — Tm)/Tm as discussed in the original 
PFC model (H. 



2. Simplification of scaled binary model 



The rescaled PFC dynamic equations can be further simplified to a lower order form via a scale analysis. 

A simple scale analysis of Eqs. ([55)) and (|34l) yields tt,, 7/1 ~ 0(6-^^^) (e.g., from Eq. ([55)) we have 0{en) ~ ©(tt,'^) and 
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0{ip) ~ 0{n), as is usually assumed). To simplify the results the following approximations are made: (i) Assume 

that {|do|, \Co\, \SC^o^^\} > {I^Cnl, jcf \ lAC^^^]} and \pf - pf\ < \pf + pf\. (An example case would be that 
the zeroth-mode {q = 0) correlation functions between the same atomic species are of the same order, and are either 
significantly larger or significantly smaller than those between different ones; see Eq. (f23|) .) Thus for the rescaled 
parameters in Eq. (f35|) . we can estimate (based on the definitions in Eqs. ([22|) and ([23|) . as well as Eqs. ([32]) and 
([M)) ) that 

go,9,g2,Vi,UQ - 0(1) or 0(e^/^), 5, 5i, 52, "i, ^2, f2, ffs, ^o, 7i, 72 ~ 0{e) or higher. 

(ii) The concentration field ip is slowly varying in space, and we can keep only the lowest linear gradient terms for ip. 

(iii) Similar to the single-component case in Sec. Ill Al it can be argued that in Eqs. ((33)) and (p4l) . compared to the 
first terms V^{ - • •}, all other terms (<7o^ ■ {■ ■ '}) of higher order, (iv) For linear terms in n, only [— e-l- (V^ -|-(7g)^]n 
is kept which will lead to periodic crystal structure in solid phases, while the other term (q;2V^ + ^ V'')n is neglected, 
which corresponds to i gno ring the nip related terms in the free energy functional (j2ip owing to to the much larger 
length scales of ip field [2a. (v) It is assumed that a2 — a4/2 ~ 2a (see the discussions below Eq. ([24]) V 

To lowest order in ©(e^'^) the above simplifications reduce the PFC equations ([3T|) . p3|) . and ([34]) to 

dn/dt = Vi+mV2, dip/dt = mVi +V2, (38) 

where 

Pi = V2 I (-e + gip^) n + (V2 + q^fn + 32"^ + + vi^J^ + 2ao [ip (V^ +V^)n+ (V^ + V'') (ni/j)] | , 

V2 = [2aon (V^ + y^)n+ (wq + 2vin + gn"^) ip + uqiJj^ - KoV'^ip] , (39) 

with ao = goCi the rescaled solute expansion coefficient. Equations (|38p and (|39l) recover the original binary PFC 
model with conserved dynamics for both n and ip fields [25l l32j . except for the vi terms {vitp"^ and 2virLip), which 
account for additional coupling between the atomic density and concentration fields (or between small crystalline and 
"slow" concentration scales). This can also be seen via rewriting Eq. (15^1) through an effective potential (or free 
energy functional) J-"eff: 

2?i=V2^, I?2 = V2^, (40) 
J-^ff = / dr j-ien^ + \n (V^ + qlf n + \g2V? + + 2ao" (V^ + V*) (nil)) 

+i(u;o + 2«in + gn2)V/ + \uqiP^ + |VV^|'} . (41) 

In the rest of this work, all results, including the corresponding amplitude equation formalism, noise dynamics, and 
the related applications, are based on the simplified PFC dynamic equations ([55)) and ([5^ . 

The above results can also be derived and verified through two other alternative methods, as given in Appendix 
\F\ Furthermore, to include higher-order terms (e.g., n^, V'^j ••■) in both the free energy functional and the dynamic 
equations, we need to consider higher-order direct correlation functions (4-point, 5-point, etc.) as shown in the 
single-component case (Sec. Ill A[) . with similar derivation steps. 



III. AMPLITUDE EQUATION FORMALISM FOR BINARY PFC MODEL 

As discussed in Sec. U the PFC methodology includes model equations governing the dynamics of density and 
concentration fields as given above. This section will examine the long wavelength and time limits of the alloy 
PFC model by deriving its corresponding amplitude equations, which emerge after coarse-graining the model using 
a multiple-scale analysis. The amplitude representation for single-component PFC models has been well established 
po'-'so'l , while for binary systems the corresponding amplitude equations have been derived very recently, for both 2D 
hexagonal/ triangular and 3D bcc and fee crystalline structures |43||. Here we focus on the 2D amplitude equations 
for the binary PFC model with hexagonal lattice structure, yielding a complete formulation incorporating the effects 
of different mobilities between alloy components and dynamic variation of the average atomic density, both of which 
are missing in the previous binary alloy amplitude formulation [32]. It is straightforward to extend this calculation 
to 3D bcc or fee structures. The derivation process involves two steps: the standard multiple scale expansion Q 
for lowest order amplitude equations (Sec. IIII Al) . and a new hybrid approach (combining results of multiple scale 
approach and the idea of "Quick and Dirty" renormalization-group (RG) method developed by Goldenfeld ei al. 
[IE 1131) for fiiU order amplitude equations (see Sec. IIII B[) . To apply the multiple scale analysis, the rescaled PFC 
equations (1551) and (15^)) are used. 
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A. Multiple scale expansion: Lowest order amplitude equations 

Following the standard procedure of multiple scale approach 85], in the limit of small e (i.e., high temperature) we 
can separate "slow" spatial and temporal scales {X = e^^^x,Y — e^^'^y,T — et) for structural profile/amplitudes from 
"fast" scales of the underlying crystalline lattice. Substituting 



and the expansions 



n = el/2„(l/2) + + g3/2„(3/2) ^ ^2^(2) ^ . . . ^ 

^ = £1/2^(1/2) + £^(1) + £3/2^(3/2) ^ ^2^(2) ^ . . . ^ 



(42) 



(43) 



into the PFC equations (|38p and (|39|) . we can obtain the corresponding equations at each order of e^/^. For sim- 
plicity, assume m,ao,g,uo,Ko ^ 0(1), 92, vi ^ ©(e^/^), and wo ~ 0(e) (as also assumed in Sec. Ill B 21 for model 
simplification). To ©(e^/^) and C(e) we have 



0, 



0, 



(44) 



where i = 1/2 or 1, and £0 = (V^ +(?o)^- This leads to (1 - TO2)y2£^^(i) ^ q ^^^^ (1 - m2)V'*'0(*) = 0, with solutions 



— n, 



e"'^ ' +C.C., 



(45) 



i=i 



where represent the three reciprocal lattice vectors for 2D hexagonal/triangular structure: = —qo{\/3x/2 + y /2), 

Q2 = 90^7 93 = qoiVSx/2 — y/2). Aj are the slowly varying complex amplitudes of the modes q^, while hq and V'o 
refer to the real amplitudes of the zero wavenumber neutral mode as a result of order parameter conservation [86j . 
Expanding to ©(e^/^) yields (with V, = (Sx, 9y), V • - d^dx + dydy, and = a| + 9f-) 



- (2V • V,)" - q^W^ 



,(1/2) 



32^(1/2)^ +„(l/2)^ 



gV^'/') ^2^(1/2) + 2aoV2 [^-(1/2) (2V • V,) n^/^) + (2V • V,) (v^^^'^r 
TO |2aoV^ 

TO I 



(1/2) 



,(1/2) 



{2V-Vs)n 



V' (2V • V,)^ - q^V 



- #(i/2)v2„(i/2)2 _ 2t;iV/^/')v2n(i/2)} , 

j(l/2) _ v2 



gV'(i/2)2y2n(i/2) + 2aoV2 ^(1/2) (2V . V,) ; 



.92" 

,(1/2) 



(1/2)2 (1/2)3 



(2V- V,) ^ 



',(l/2)„(l/2) 



)]} 



9^^,(1/2) + 2aoV^ 



,(1/2) 



(2V- V,)n 



(1/2) 



g^(i/2)v2„(i/2)^ _ 2i;i^(i/2)v2n(i/2), 



which is equivalent to 
(1 - TO2)v2£on(3/2) = 



drn^'/^^ ~ mdr^P^'/^^ + (1 - m^) | [v^ - (2V • V,)' - g^y^j n(V2) _ 5^(1/2)^2^(1/2) 

■^(1/2) (2V . V,,) + (2V • V,) } , 



■g2n(i/2)%„(i/2)3- 



2aoV 



(1 - to2)KoVV^'/'^ = TO^Tnti/^) _ aj,^(i/2) 



(l-m2){- 



2aoV 



,(1/2) 



(2V- V,)n 



(1/2) 



■ gV^(i/2) V^nd/^)' + 2i;i^(i/2)^2^(i/2) I (4g) 



As shown in Eq. (|l5|) . the zero eigenvectors of operators V^/Iq and are (e^"'3 '',l) and 1 (of the 0th mode), 
respectively. Using the Fredholm theory or solvability condition i85;] in Eq. (|46p . we can derive the lowest order 
amplitude equations as (with j = 1, 2, 3) 



-1 + (2^gO • V, 



-3A 



(1/2) 



A 



(1/2) 



2E 



fc<i 



A' 



(1/2) 



A 



A 



(1/2) 



(1/2) 



o (1/2)2 (^y2) , ,(1/2)^ 



(1/2) 
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-2ao [4^/^^ {2^g^ ■ V.) ^^'^ + {2^q^ ■ V.) {i^'i"^ A^'^)] } , (47) 

dn^^"^/dt = qtVln^^"\ (48) 
d^''^^'^^ /dt^mq^Vl^'^^'^\ (49) 

Using the scaling relation Aj = e^/'^ A'"^^'^\ uq ~ e^/'^n'^^^\ and tpQ — e^^'^'ipo^^'^\ we can then obtain the corresponding 
amplitude equations in the unsealed units (x, y, t). 

It is noted that the direct solutions to Eq. have the form 

3 

„(3/2) ^ ^(3/2) Y,T)+J2 Af^^'^ (X, Y, T)e*95'-^ + c.c. + higher harmonics, (50) 

3 

^p(3m ^ ^jf^^ [X, Y.T)^Y. Y, T)e"iy'' + c.c. + higher harmonics. (51) 

Compared to Eq. for the 0{e^/'^) and 0{t) solutions, it can be found that the complex amplitudes ipj corre- 

sponding to the periodic modes of the concentration field in substitutional binary alloys considered here is generally 
of order e higher than Aj^ uq, and tpQ. For systems in which a sublattice ordering occurs (such as B2 or B32 ordering 
in bcc crystals), V'j would be of the same order as these other fields. To describe sublattice ordering a different free 
energy functional from the one given in Eq. (1411) would also be required. Detailed results will be presented elsewhere. 



B. A hybrid approach: Full order amphtude equations 



The lowest order amplitude equations p7l) - p9)) derived above are not sufficient to describe the evolution of binary 
systems; e.g., Eq. (H7)) for Aj is not rotationally invariant, and Eqs. (051) and (05]) for no and ipo are just diffusion 
equations and would lead to a steady state solution of constant no and V'o values at long enough time. We thus 
need higher order amplitude equations, which in principle can be derived by extending the multiple scale process 
described above to higher order expansions. However, the procedure is complicated and tedious. In the following we 
use, instead, a simplified appro ach combining the above steps of multiple scale expansion and the idea of the "Quick 
and Dirty" RG method [11, l^. 

The first step is the standard multiple scale expansion given in Sec. IIII Al starting from the scale separation Eq. 
(021) ■ From the solution forms of Eqs. (HSl) . ([501 and (ISTI) . we know that to all orders of e the solutions of n and ip 
fields can be written as 



n = no{X, Y,T)+J2 Aj {X, Y, T)e'''^ + c.c. + higher harmonics, 

3 . 

^ = 'ipo{X,Y,T)+ J2 ipjiX,Y,T)e"ir'- + c.c. + higher harmonics. 



(52) 
(53) 



with {X,Y,T) the slow scales. Thus, based on the separation between "fast" /" slow" scales the following expansions 
(full-order) can be obtained: 



V^n -> eV2no + ^(/:j^Aj)e*'^?"' + c.c. + {•••}, 
{V' + ql)'n ^ {eWj + ql)\o + Y.iSfA,)e"'"- + c.c. + {■■■}, 



(V2 + V'')(ni/;) -> {e^l + e^\/i) ino^po + J2'^i'^3 



E 



k<l 



e"ir- + c.c. + {■■■}, 



k<l 



-> ng + 2^ |^,f + ^ 2no^, + 2 ^ AlA^ | e^^^r + c.c. + {■ ■ •}, 
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n ^ n 



g + 6no^|A,f + 6 Yl A, +C.C. 

3 ( k<l 

+ { 3("o + 1^. + 6 E M^^r + a, + l e^<^ + c.c. + {. . .}, 



fc<i 



2no ^0^, + ^*k^i 



+ 2 ^ iV'fcp A, + 2^0 ^ A* vr + 2v^, Y (Aurk + c.c.) + a*^ 



■'■ + C.C. + {•••}, 



(54) 



where {■ • •} refers to the contributions from higher harmonics and the slow operators are given by 

q = eVl + (2,qO . _ ^2^ gs ^ ^ ^2 ^ ^^2 ^ ^1/2 (2,gO . (55) 

Assuming that higher harmonic terms can be neglected, the binary PFC equations (p8l) and ([39]) are then replaced 

by 

edrno + SrAj-e^'j" '- + c.c. = Pf + ml?|, (56) 

eSrVo + e Ej SrV'je'''^"-'' + c.c. = mPf + 2?|, (57) 

where 2?f and I?| are the corresponding expansion of I?i and I?2, as obtained by substituting Eq. (|54l) into Eq. ([39|). 
Integrating Eqs. ([55]) and ([57)) over the eigenmodes / drje^^'^J 1}, keeping in mind that "fast" and "slow" scales 
are separated, and in the final step returning to original unsealed units {x,y,t), we arrive at the following full-order 
amplitude equations for the binary PFC model: 



dno dt = V^— + mV2 — . 



dAj/dt = Cj 



5F 5F 



o 5F ST 



^ , ,^ ^ f 6T St\ 2 f ST ST 



(58) 
(59) 
(60) 
(61) 



where j — 1, 2, 3, and 

T=\dT 



{11 11 ^ 

3 3 / 3 \ 

+ E l^.^.f + I E l^^l' + + 252) n + c.c. + 6^ |A,f 



+5 



3 



j<fe 



^"oV^o + "0 E l^^l' + 2 E + E ( 2"oV'oA,^* + \A]i,f + c.c 



j,fc=i 



^^(A,V*+c.c.)(A-^:. 



+2ao 



- c.c.) + Y ("o^j + V^o^j )^fcV'r + c.c. 
3 



i>ono (V^ + V'') no + V-o I E ^^t^fijA^ + c.c. 
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+nQ^il;*CjgjAj + ^ AjtpkCiGiAi + c.c. 



3 3 

-lyoV'o + T^Ko |V?Ao|' + tWoV4 + {wo + Suo^o) ^ I'Ajf " ^^^o ^ (^j-CIV* + ex. 



+ Vl 



j<k 



Corresponding to Eq. (I55|) . the operators £j and CJj (in the original scales) are defined by 



(62) 



(63) 



and for simplicity, in Eqs. (j59p - (|62p the operator Lj can be replaced by —q^ in the long wavelength approximation 
as adopted in Ref. j3^ . 

As discussed at the end of Sec. IIII A[ the amplitudes i/'j are of 0(e) higher compared to the others for the free 
energy functional considered here. Thus the above amplitude equations can be further simplified by assuming 7/;^ ^ 0, 
which leads to 



dAjjdt 



ST 



mql I 2aQ 



A, (V^ + V-*) no + no£,g,A, + ^ Al^G^A^ 



k<l 



2gVo(no^j + AlA*i) + 2viijoAj 



ST 



ST 



The dynamic equations for uq and ipQ are still governed by Eqs. (j58p and (j60l) . The amplitude equations can be 
further simplified by noting from Eq. (|5T|) ~ dij^j/dt = —qQ{'mST/SA* + ST/S'ip*\^.^o). Thus, the above dynamic 



equation for Aj can be further approximated as 



dAj/dtc^ -qlil-m^) 



ST 

SA*' 



(64) 



which to lowest order recovers the result of multiple scale approach given in Eq. (|47p . In the applications that will 
be examined in Sec. |V]the simplified amplitude equations ([58|. (|60| . and (|64| will be used. 



IV. NOISE DYNAMICS AND STOCHASTIC AMPLITUDE EQUATIONS 

In the original PFC model [1^, [13] a conserved noise dynamics has been incorporated. However, in DDFT it has 
been argued that the dynamic equation governing the density field evolution should be deterministic and an additional 
stochastic noise term added to Eq. ([T]) would lead to an artificial double-counting of thermal fluctuations f76||. On 
the other hand, recent studies [87j have shown that deterministic DDFT dynamics governs the ensemble averaged 
density field p{r, t), while if the density field is temporally coarse-grained -as is the assumption in PFC modeling- the 
corresponding dynamic equation would then be stochastic, but with a (unknown) coarse-grained free energy functional 
instead of the equilibrium Helniholtz free energy functional used in static DFT. In the current case of PFC modeling, 
quite drastic approximations have been made to the DFT free energy functional (particularly at the level of the 
direct correlation functions; see e.g., Eqs. ([3]), (Hj), and ((TSl) ). and hence it could be argued that the incorporation of 
noise terms in the PFC dynamics would be necessary and useful to capture the qualitative effects of fiuctuations in 
phenomena such a homogeneous nucleation. In what follows, noise will be added to the PFC models studied above 
and the corresponding stochastic amplitude equations will be derived for both single component and binary systems. 
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A. Single-component PFC 



The stochastic DDFT equation for single-component systems is given by Eq. ([T]) with a muhiphcative noise term 
V • [\/ pit", t)C{r, t)] added to the right-hand-side, where the noise field C('"i^) is determined by (with Fq = 2kBTM) 

(C(r,t))-0, {C>^{r,t)C{r',t'))=T^5{r-r')6{t-t')5^'^ {fi,,. ^ x,y, z). (65) 

The corresponding dynamic equation governing the rescaled density field n is similar to Eq. ([7]), i.e., dn/dt = 
M'V • [(1 -I- n)'VSJ'/Sn] + V • [V (1 + n)/pi Adopting the lowest order approximation as given in Sec. Ill A| we 
can write the rescaled stochastic PFC equation as 



dn/dt = [-en + (V^ + qlfn - 



.92"- 



+ V • c, 



(66) 



where the rescaled noise C is also determined by Eq. (|65|) but with Fq ~ 2v / [B^"^ pi) (where d is the dimensionality). 

To derive the associated stochastic amplitude equations, we follow the standard multiple scale approach in the limit 
of small e, which leads to the expansion of density field n in terms of the zeroth-mode average density ng and complex 
amplitudes Aj that are varying on slow scales {X,Y,T)\ see Eg. ([52l). Effects of external noise can be approximated 
via a projection procedure used in hydrodynamic analysis (ssl . l89l| . Based on the fact that thermal noises originate 
from the fluctuations or random motion of individual atoms/molecules at the microscopic scales, we can project C, 
onto the base modes given in Eq. ([5^ . i.e.. 



C = Co(^, Y,T) + Y, Ca, {X, Y, T)e<-^ + c.c. 



(67) 



where 



(Co) = {Ca,) = 0, {CaXa,) - (CoCa,) = (Coa,> = 0, 

(Co^Co") = d^T^Sir - r')5{t - t')5^\ (ClG,*) - ^.To5{r - r')5{t - t')5,,5'^^ 



(68) 



(with j, J = 1, 2, p,v — X, y). Here iJj (z = 0, 1, 2, 3) is a constant determining the noise correlation strength, which 
can be approximated as di = d — 1/1 \i equal contribution from all modes in Eq. (|67p is assumed. Thus the random 
noise term in Eq. (j66p is given by 



1/2 



(69) 

©(£3/2) and 



In order to be relevant in the amplitude expansion, it is necessary that V • C ^ 0{e^/'^), leading to C,Aj 
hence the noise intensity Fq ~ C(e)- The latter yields ^ C'(e'^/^), which can be deduced from Eq. 

Following the procedure of multiple scale expansion and retaining the random force contribution to the lowest order, 
we can derive the following stochastic amplitude equations 

dAj/dt = -ql6F/6A* + Q , (70) 

dno/dt = W^ST/Sno + V • Co, (71) 

where is the effective free energy of the single-component amplitude representation (sec Refs. [1, [s^, [13] for the 
detailed form), which is given by Eq. ([5^ with ■00 and ipj set to 0. Also, (j = iq^ ■ Caj (j = 1, 2,3) and 

(0) = (Co) = 0, (GO) - (CoO) - (CoC;) = 0, 



(GC;) = ^^q^o^o5{r ~ r')5{t ~ t')%, (Co^Co") = ^9oro<5(r - r')8(t - t')5^^- 



(72) 



The noise dynamics is then consistent with the dynamics of amplitude representation, i.e., nonconserved dynamics 
for Aj in Eq. (|7D|) and conserved one for no in Eq. ([7T|) . 



B. Binary PFC 



Similar to the single-component system, based on Eq. (jlSp the stochastic DDFT equations for a binary system can 
be written as 



dpA 
dt 



V ■ 



MapaV— + vpICa 

OPA 



dpE 
dt 



= V 



MbPb"^-. h ^/pbCb 

dps 



(73) 



15 



where for noises (a, P = A, B; fi,v = x, y, z\ Y a = 2kBTMa)^ 

(0(r,t)) =0, {(:^^{r,t)Cp{r\t'))^T^5{r-r')5{t-t')5^p5'^r 
From Eqs. (jl9p and (|26p the dynamics equations for n and fields can be rewritten as 

dn/dt : 



(74) 



MiVi + M2V2 
1 



V 



ViT^ v/i + ^Ca + Vi^Ci 



1 



{{M2 - AfiV)2?i + (Ml - M2V')2?2 



(75) 



+(1 - V)V • [v/(l + n)(l+^)CA] - (1 + V^)V • [V(l+«)(l-V')Cs] } 



where we have rescaled Cyi(B) ~^ Ca{B)/V^Pi- Following the procedure given in Sec. Ill Bl and only retaining the lowest 
order noise terms, we can derive the rescaled stochastic binary PFC equations as 

dn/dt ^Vi + mV2 + V ■ Cn, dip/dt ^ mVi + V2 + V ■ d,, (76) 

where the expressions of 2?i and 2?2 have been given in Eqs. (I5ni) - (PT|) . The noise terms are defined by 

Cri = Ca + Cb, Ci> = Ca - Cb, (77) 

where Ca{b) also obeys Eq. (|74|) . although with Fq — ksT MaV / {MiBq"^ R'^ pi) due to the rescaling, and 

(Cn) = (a) = 0, (C^C;;) = (r^ - TB)5{r - r')S{t ~ t')6^\ 

{QO = (C^Cp = (r^ + rB)S{r ~ r')S{t ~ t')S^-^, (78) 

with F^ + Fb = 2v/{BfR'^pi) and F^ - F^ = m{TA + F^) = 2mv / {Bf R'' pi) . 

Using the multiple scale approach, we can expand the density field n according to Eq. (j52p while assuming the 
concentration field as slowly varying, ip — tpQ(X,Y,T) (that is, keeping only the zeroth mode and neglecting the 
higher-order contributions from ipj in Eq. (j53p . as discussed in Sec. IIII Bp . Similar to the single-component case, the 
projection of noises can be given by 



C„ = Co(^, Y,T) + Y, Ca, {X, y, r)e*'J?- + C.C., Cv = U,{X, Y, T). 

Thus the expression of V • Cn is the same as Eq. (p^ . while V ■ Cip = ^^^^{dxC^ + dyC^) 

Ca, , Co, Q ~ 0(e'/') and F^, Fb ~ 0{e). 

The stochastic amplitude equations for binary PFC model can then be derived, i.e.. 



dno/dt 



Sno 



ST 



dno 



, ST 

Sijjo 
, ST 

SlJjQ 



V-Co, 



(79) 

Also we can estimate 

(80) 

(81) 
(82) 



where the deterministic parts have been obtained in Sec. IIIIBj see Eqs. ([58| . ([60| . and (|64p . as well as Eq. (|62|) for 
the effective potential T. For the noise terms, Q = iq'j ■ Ca, {j — 1, 2, 3), and 



(0) = (Co) = {QJ = 0, (CO) = (CoO) = (Coc;) = (c^oO) = (c^oC;) = o, 

(C.C;) = ^^qU^A + rB)S{r - r')S{t - t')%, (Co" Co") = M^a + rB)(5(r - r')S{t - t')S^- , 

(CC^o) = (r^ + rB)S{r - r')S{t ^')'5''^ (C^„Co") = (r^ - rB)S{r - r')Sit - t')S^'^, 



(83) 



with i,j = 1, 2, 3 and /i, i/ = x, y. If assuming Ta — F^ (for equal mobility AIa — Mb and m ~ 0), i.e., with almost 
the same noise/fiuctuation intensity for A and B components, we have (C^qCo) — ^ and hence all noise terms (Cj, 
Co, Cifia) can be treated independently. However, for the case of different mobilities {AIa Mb and m ^ 0), we get 
iCipgCo) 7^ 0, and hence noises Co and Cv'o for 0th- mode density fields uq and ipo are then correlated. Similar results 
can be obtained for noises Cn and Ci/i in the stochastic PFC equations ([TS)) and ([75]) . 
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V. APPLICATIONS IN ALLOY HETEROSTRUCTURES 



As discussed in the introduction, the PFC model and the corresponding amphtude equations have appUed to 
the study of a wide variety of phenomena involved in material processing and microstructure evolution. In this 
section we will illustrate how the amplitude equations derived in the preceding sections can be employed to examine 
the effect of surface segregation and alloy intermixing. Alloy intermixing is known to play an important role in 
the growth and processing of material heterostructures, including morphological and compositional profiles and the 
associated sample optoelectronic properties and functionality. Recent intensive studies on thin film epitaxy and atomic 
deposition have shown the important effects of intermixing on nanostructure self-assembly. Typical examples include 
InAs(InGaAs)/GaAs(001) or Ge(SiGe)/Si(001) heteroepitaxy that has been investigated extensively 

(particularly the intermixing-caused alloying of wetting layers and quantum dots), and the interlayer diffusion in 
semiconductor multilayers or superlattices such as InP/InGaAs [7l|, GaAs/GaSb [zl], and GaAs/InAs [73]. An 
important phenomenon in these epitaxial layers is the occurrence of surface segregation, in which an enrichment of 
one of the film components at a surface or interface region occurs. This has been observed in a variety of material 
systems including III-V and II- VI semiconductor heterostructures [66l - [73j . To address these complicated phenomena 
and effects, the basic processes and mechanisms of intra- and inter-layer diffusion at nearly-planar interfaces as well 
as their coupling with material processing and growth parameters needs to be clarified. 

In light of the above observations, the focus of this section is on heterostructures of a nearly-planar interface, for 
both lattice-matched and strained epitaxial layers. For layers stressed due to lattice mismatch, the configurations 
studied here are metastable in nature, and our results will be used for further studies of the associated later-stage 
nanostructure evolution (e.g., quantum dots), which will be presented elsewhere. For such film geometry with a 
planar interface, it can be assumed that both morphological and compositional profiles along the lateral direction 
are approximately uniform or homogeneous (at least metastably), and hence these structural profiles vary only along 
the direction (y) normal to the interface. An advantage of the amplitude equation representation of the PFC model 
developed above is that the system of interest can then be mapped onto an effective one-dimensional (ID) description, 
as will be shown below. 



A. Effective ID model system with elasticity 



To address the elasticity incorporated in the amplitude equation formalism, it is useful to note that the structural 
amplitudes can be written as 



A, 



A'e^"." (j-- 1,2,3), 



(84) 



where for 2D hexagonal structure are the three basic wave vectors given in Sec. IIIII and u — Sq(xx + yy) describes 
the bulk compression or dilation. The effective free energy J- in Eq. (|62p can be rewritten as (neglecting the higher 
order contributions from ipj and approximating Cj ~ — ^q) 



where 



E + 1 E 141' + (6«o + 2.92) 114 + c-c- + 6E I4n4i 



+2ao 



V^ono (V^ + V*) no - 9oV^o E 4*^j4 + ^-C- 



2i 6, 



2g".<5, 



(85) 



(86) 



with Si — —SxX — Syy/2, 82 = Syij, S3 = SxX — Syy/2, S^ = a/39o'5o/2, and Sy = qoSq. The corresponding dynamic 
equations for A'^, no, and ■00 are still governed by Eqs. (|5n |) - ((5^ . although with Aj replaced by A'j. In mechanical 
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equilibrium, we can assume that A' ~ A, i.e., Aj 



Aexp{iq'^ ■ u) where A is a constant. Minimizing the effective 



free energy IF with respect to A yields the equilibrium value (5q 



-1 + y/1 — 2aoV'o — — aoV'o to lowest order. This 



leads to the equilibrium wave number goq = (1 + '^o'')'Zo = Vl ~ 2q;oV'o 9o (where ao is the rescaled solute expansion 
coefficient defined in Sec. Ill Bp . and the equilibrium amplitude 



A 



^ |-(3?io + 92) + \l (3no + <72)2 - 15 [-e + q^{5l + 2Jo)(5g + 25^ + 4aoV'o) + no(3no + 252) + <?V'o] 



(87) 



The elastic constants (rescaled) are then given by Cn = C22 = 9A^, C12 = C44 = Cn/S = 3A^, and Young's modulus 
E = 8^2 [13, [11,133. 

For the dynamics of a heterostructure configuration with nearly-planar interface (either liquid-solid or solid-solid), 
we can assume that A'^{x,y,t) ~ A^j(y,t), nQ{x,y,t) ~ ng(y,t), and ipo{x,y,t) ~ ipQ{y,t), resulting in an effective ID 
description of the system. The dynamics of the amplitude equations then become 



dn°/dt 



dA^/dt = -qlil-rn')^, 



(89) 
(90) 



where 



6T 



6^ 

Hi 



+2ao [Vo° {d^ 4 
(wo - Kodl)4,° + u^^f + g \ 



df, I ri, 



nl + g^nf nf + (Gn^ + 252) ^ + 6 | H ^" + 





y) "-0 



'2 1 ^0 + 2z;in°Vo° 



(91) 



5T 



-2ao 



-e + gf + 2g2nl + inf + g^jf A° + 



-(6nSJ + 252) n At - 2a,ql K^^A" + g° (^{JA^)] 



(92) 



1^012 + 2^(1^2^2 + 1^?^ 



(93) 



with 



2qlS,. 



(94) 



For coherent strained alloy layers, which are of great interest in materials growth, the solid layer is strained with 



respect to a substrate and subjected to an epitaxial condition 



„sub 



(\/3/2)(jo(l + (5g"^) (with "sub" referring 



to the substrate). The wavenumber qy along the vertical or layer growth direction y is determined by the lattice 
elastic relaxation (or Poisson relaxation in continuum elasticity theory). The system is thus governed by the above 
amplitude equations (|88l) - ([93|) . but with Sq fixed by the corresponding elasticity quantity (5™^ of the substrate (and 
thus Sx — ^/SqaSfl^^ /2 and 6y — go^o"'')- The vertical strain relaxation (Poisson relaxation) can be determined from 
the phase of complex amplitudes A^. Furthermore, the misfit strain of such a solid layer is given by 



R 



R 



1 



So 



"0 



(95) 



where R and q^ are lateral lattice spacing and wavenumber of the strained layer, and Rcq, qx,cqi and ^q'^ are for the 
corresponding stress-free, equilibrium bulk state. 
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For the systems studied here the model parameters are chosen such that no phase separation or spinodal decom- 
position can occur in the bulk of each solid or liquid region. The corresponding conditions on the parameters that 
assure this are derived via a linear stability analysis of the amplitude equations. Following standard procedures, we 
substitute the expansion rig = + '^Oi "00 = V'o + V'Oi and Aj = Aj + Aj into Eqs. (155)) . (|Bn|) and ([M]) . obtain the 
linearized evolution equations for the perturbed quantities no, V'o, and Aj, and calculate the associated perturbation 
growth rates. The corresponding results are complicated due to the coupling between the evolution equations of all 
three perturbed quantities. To estimate the conditions for phase separation, here we simply assume that no,Aj ^ 0, 
and only study the stability of concentration field. To first order of V'o we have 

dipo/dt ~ j-isToV^ + wq + SmoV^o + 9^1 + 2wino + 2gY^\Ajf + m [2aono(V^ + V) + 2gno-0o + 2ui-0o] | V'o- 

(96) 

In Fourier space, the perturbation growth rate a{q) is then given by 

a = -q"^ [2maonQq'^ + {Kq - 2maQnQ)q'^ + w^e] , (97) 

where 

Wcs = wq + Swo'/'o +5^0 + 2vino + + 2m{gno + wi)V'o- (98) 

j 

If WcB < 0, an instability of the homogeneous alloy occurs, leading to spinodal decomposition or phase separation 
of alloy components. The characteristic wave number (for maximum perturbation growth rate) is then given by 
(Jmax = W (Ko - 2maQno)^ - GmaofioWcS ~ (Kq - 2maono)] / {Gmaono) if m,ao,no ^ 0, or ql-^^^ = -Wcs/i2Ka) if 
one of m, "-o = 0. 

For the heterostructural systems presented below and the parameters chosen, the condition Wctf > is always 
satisfied in the bulk phases, keeping homogeneous concentration profile within each layer. Concentration heterogeneity 
may occur across the system configuration, which however is due to the effect of interfaces or due to composition 
overshooting, a phenomenon caused by alloy intermixing that will be discussed in detailed below. 



B. Results: Equilibrium profiles and layer growth 



Equations ()88p -(|93 |) were solved numerically using a pseudospectral method and an exponential propagation scheme 
for time integration of stiff equations [90,[9l| ■ Results of the corresponding morphological and compositional ID profiles 
are shown in Figs. [iHll for two types of configurations of liquid-solid-solid and liquid-solid coexistence or growth. For 
the simulations shown here we choose a time step At = 1, which can be made as large as this due to the numerical 
scheme we used; The numerical grid spacing used is Ay — Xq/S (where Aq = 27r/qo)- To emulate a liquid-solid (or 
liquid-solid-solid) heterostructure and apply periodic boundary conditions in the numerical calculation, the initial 
configuration is set as two (or four) symmetric interfaces located at y = Ly/A and iLy/A (or y — Ly/Q, Ly/i, 2Ly/?> 
and 5Ly/6), separating different liquid or solid regions. These interfaces need to be set sufficiently far apart from 
each other to avoid any interface coupling and the artifacts of finite size effects. For results shown below we choose 
the ID system size perpendicular to the interfaces as Ly = 2048Ay, with similar results obtained in calculations up 
to Ly = 8192Aj/. Also, the parameters used in the amplitude equations are based on the phase diagrams given in 
Ref. [s^l showing liquid-sohd and sohd-solid coexistence, i.e., (g, (72, "o, ^o, wi) = (—1.8,-0.6,4,1,0), wo = 0.008 or 
0.088, ao = 0.3 or 0, and e = ±0.02. 



1. Liquid-solid and liquid-solid-solid coexistence 



The equilibrium profile for a liquid-solid(I)-solid(II) coexistence is given in Fig. [T] (with time corresponding to 
t — 2x 10''). To obtain the liquid-solid-solid coexistence, we use e = 0.02, ao = 0.3, and wq = 0.008 (from the eutectic 
phase diagram in Ref. [l^l), set the initial length ratio of liquid:solid(I):solid(II) as 1/3:1/3:1/3, and let all of ■00: 
and Uq evolve with time until a stationary state is reached. Solid II is treated as a substrate (unstrained), and hence 
in the amplitude equations (j88| -(|93 l) we set Sq = S}^ = —1 + ■\/l — 2aoV'o^- Due to nonzero solute expansion coefficient 
ao, i.e., different atomic sizes of A and B alloy components, solid I is strained (with misfit £,„ with respect to the 
substrate (solid II) being 14.9% for the parameters of Fig. [T]). This is consistent with the numerical results in Fig. [TId, 
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y/Ay (perpendicular to the interface) y/Ay (perpendicular to the interface) 

FIG. 1: Liquid-solid-solid coexistence profile calculated from the amplitude equations, as characterized by (a) the composition 
field ipo, amplitudes \Aj\, and the average density field no, and (b) the phases of amplitudes Aj. The parameters are set as 
e = 0.02, Qo = 0.3, {g, g2,uo, KojWo) = ( — 1.8, —0.6,4, 1,0.008), and the time at t — 2 x 10^. Solid I is strained with respect to 
the substrate solid II. 

showing zero phase of ampUtudes Aj within unstrained solid II and a linear dependence of phase on position y in the 
bulk of solid I. For comparison, the magnitude of lattice misfit between III-V or II- VI layers is around to 5% (e.g., 
Em = 4.2% for Ge/Si and less for SiajGei-^/Sij^Gei-j,), while the lattice mismatch for III-V Nitride heteroepitaxial 
films or Ill-V/Si heterostructures could reach 10% or more (e.g., Em = 11.5% for InAs/Si). 

For a liquid-solid heterostructural configuration, to determine the coexistence state we choose similar parameters 
except for wq = 0.088, e = —0.02, and initially ^/jq = in the whole system. This corresponds to the single solid 
phase region (no solid-solid coexistence, only liquid-solid) in the phase diagram 32]. To make the solid strained, we 
set Sq = 0.05 as given by an external condition (i.e., a substrate), and thus from Eq. (|95p the misfit strain in the solid 
here is about 5%. The results for uq — and 0.3 are given in Figs. [5^ and[2)D respectively, including the equilibrium 
profiles (up to t — 2 x 10^) and the process of time evolution. As expected, for ao — (equal atomic size of alloy 
components) the concentration field ipo remains at all the time, as seen in Fig. [5^. However, for ao = 0.3 the 
initial V'o = profile splits at the liquid-solid interface (see Fig. [^b)- For the parameters used here, ao > (with 
size of atom A larger than that of atom B) and misfit > (compressed solid), and thus the solid would prefer to 
have more smaller atoms B (with -00 < 0)i leading to a "dip" on the solid side of the compositional interface; due 
to the conservation law on the field V'Oi a "bump" of "00 > (more larger atoms A) appears on the other side via 
layer intcrdiffusion or alloy intermixing. As a result of atomic diffusion, such "dip" and "bump" will spread out into 
the bulk phases as time increases, leading to a positive/negative 00 equilibrium profile of liquid-solid coexistence, as 
shown in Fig. [^h- 

It is interesting to note that the non-homogeneous compositional profile can also be found in liquid-solid heterostruc- 
tures with nonzero ao and no misfit strain (i.e., 6o = 0, ipo = 0, and ao = 0.3, as in Fig. [3]). A slight enrichment 
of larger atoms A is observed on the surface of unstrained solid, showing as a "peak" (with ipo ~ 1.5 x 10""*) at the 
compositional interface in Fig. [31 Note that this phenomenon of weak surface segregation persists in the equilibrium 
or stationary configuration (as tested up to i = 10*"), and is caused by unequal atomic sizes of alloy components. 
Due to the conservation of the "00 field and the appearance of concentration "peak" at interface, the bulk values of 
concentration field 0^0 in both liquid and solid regions deviate from the value in the corresponding phase diagram, 
as mediated by the alloy diffusion process. We find that this deviation is a result of finite size effect: The deviation 
decreases with increasing system size, as confirmed in our simulations of Ly = 1024A?/, 2048Aj/ and 8192Ay. Thus 
in the thermodynamic limit (with Ly — > oo) V'o = is expected in the liquid and solid bulks, consistent with the 
equilibrium phase diagram for unstrained systems. On the other hand, the effect of surface enrichment would be 
preserved, as we have observed in simulations of various system sizes. 
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FIG. 2: Liquid-solid coexistence profiles, for ao = Q (a) and 0.3 (b). The other parameters are the same as Fig. [T] except for 
e = —0.02, Wo — 0.088, and initially ^0 ~ 0. The solid region has a misfit strain of around 5% due to So = 0.05 set in the 
amplitude equations. 
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FIG. 3: Liquid-solid coexistence profiles for unstrained solid layer (with So = 0) and ao = Q and 0.3. The other parameters are 
the same as Fig. [2] In the lower panel the \ Aj\ and no profiles overlap for ao = and 0.3. 



2. Coherent strained layer growth and front motion 



To simulate the process of strained layer growth encountered in most experiments, we start from a liquid- 
solid(strained) coexisting configuration and let the liquid solidify, leading to a growing front of the strained solid 
layer (as shown in Fig. 2]). The initial condition is set as the liquid-solid coexistence profiles given in Fig. [5J with 
only no in liquid changed to ng"' = —0.0021 to initialize the solidification and growth while all others (including con- 
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y/Ay (perpendicular to the interface) y/Ay (perpendicular to the interface) 

FIG. 4: Growth of strained solid layer from a liquid-solid initial configuration, with ao — (a) and 0.3 (b) and equal mobility 
Ma = Mb- The parameters are the same as the corresponding liquid-solid coexistence state given in Fig. [2l except for 
no = —0.0021 in the liquid region where a constant flux boundary condition is set up. 



centration tpo and amplitudes Aj) being kept the same as the coexistence condition. The growth rate of the strained 
layer can be controUcd by the setting of liquid Uq"^, i.e., its deviation from the equilibrium or coexistence value. A 
boundary condition of constant flux is kept in the liquid region (with distance 100 Ay beyond the moving interface). 

The growth process is shown in Fig. 31 for equal mobility Ma = Mb, 5% misfit strain for solid layer, and up to 
t — 10^. The liquid-solid front moves smoothly for both ao = and 0.3, as seen from the amplitude and no profiles 
in the figure. For ao = 0, the concentration TpQ in both liquid and solid layers remains uniform at the initial value 0, 
as in the equilibrium state. However, the results for ao = 0.3 show a phenomenon of composition overshooting at the 
growth front of strained solid (see Fig. (Ha). Such overshooting effect reveals as the increase of ipo (i.e., more A or less 
B atoms) around the interface, resulting in the phenomenon of surface enrichment: The A atoms (with larger atomic 
size for ao > 0) are segregated on the solid surface with compressive strain. As time increases, such variation of alloy 
concentration will propagate into the bulk of solid layer as a result of atomic diffusion (note that the concentration 
of liquid bulk remains unchanged due to the constant flux boundary condition) . 

Figure O shows that the mobility disparity between different alloy components plays an important role on this 
overshooting effect. Atoms with larger mobility will accumulate on the surface, even with ao = 0. As seen in the 
concentration profile of Fig. [5^, a peak of larger (or smaller) ipo appears around the liquid-solid interface for Ma > Mb 
(or Ma < Mb), while no overshooting is observed in the case of equal mobility. For nonzero ao (Fig. [SJd), the effect 
of surface enrichment of A atoms will be enhanced when Ma > Mb , while when Ma < Mb the B atom enrichment 
is observed at large enough time. 

Another effect of mobility difference presented in Fig. O is the change of solid layer growth rate or front moving 
speed. For large disparity of atomic mobility between A and B components, one of the components moves much 
slower compared to the other one and thus would hinder the atomic diffusion process. This leads to a slower motion 
of interface, as seen in Fig. [5l Thus we can expect that in the limit of Ma/ Mb ^ 1 (or Ma/ Mb <C 1), B (or 
A) atoms would be almost immobile compared to A (or B) and hence would pin the interface location, resulting 
in a frozen front. This has been incorporated in the amplitude equations developed above: When m — ±1 (with 
771 = [Ma - Mb)/{Ma + Mb) as defined in Eq. ((321)), Eq. ([64]) yields dAj/dt = 0, a frozen ampfitude profile. 
Furthermore, the concentration profile is symmetric with respect to the sign of m (i.e., Ma/Mb > 1 vs. < 1) for 
ao = 0, as shown in Fig. [5^ for Ma/Mb = 100 and lO^'^ which yield the same front moving rate and the same Aj 
and no profiles. The situation for nonzero ao (different atomic sizes) is more complicated. In our calculations of Fig. 
[5)3 with ao = 0.3 and 5% compressive misfit, the liquid-solid coexisting profile yields -00 > (A-rich) in the liquid 
region and < (B-rich) in the solid layer (see also Fig. [^b). When Ma — lOOMg, the segregation of fast A atoms 
around the interface would tend to hinder the growth of B-rich solid layer, while for Ma = Mb /WO the accumulation 
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FIG. 5: Growth of strained solid layer from a liquid-solid initial configuration, with ao = (a) and 0.3 (b), mobilities Ma = Mb, 
Ma = lOOA'/s, and Ma ~ Ms/100, and time t — 10^. Other parameters are the same as those in Fig. [l] In (a) the \Aj\ and 
no profiles for Ma = lOOMs and Ma = Mb /WO overlap. 

of fast B atoms will naturally be accompanied by the expansion of solid region, resulting in a faster solid growth. 

The composition overshooting effect presented here and the associated surface enrichment phenomenon can be 
viewed as a result of interface intermixing process via atomic interdiffusion and mass transport of alloy components, 
showing as the vertical phase separation or segregation in the liquid-solid interface region. Such process of vertical 
separation has also been found in 2D simulations of binary PFC equations [1^, where the component of greater 
size or larger mobility was found to accumulate near undulated solid surface in a liquid/substrate epitaxial system. 
Importantly, the results shown here are consistent with recent experimental observations of surface or interface segre- 
gation phenomenon in alloy heterostructures, particularly in semiconductor epitaxial layers. Most experiments focus 
on III-V or group IV heteroepitaxial films, with typical systems including InGaAs/GaAs(001) (with In enrichment or 
segregation |66.-68]), Ge(SiGe)/Si(001) (with Ge segregation [S^lzOl), and multilayers or superlattices of InP/InGaAs 
(with excess InAs at the interface ItiI ). GaAs/GaSb (with Sb segregation and Sb-As exchange and intermixing [t^). 
GaAs/InAs (with In segregation [73|), etc. In these experimental systems the segregation or enrichment effect in- 
volves the coupling of various factors of different atomic size (nonzero ao), misfit strain, and unequal mobility of alloy 
components (e.g., A/qc > Msi and Mjn > Afca), each of which has been identified in our analysis given above. 

VI. CONCLUSIONS 

In this paper we have furthered the development of the phase-field-crystal methodology by systematically deriving 
the PFC dynamic model equations from dynamical density functional theory (DDFT) and completing the derivation 
of the corresponding amplitude equation formalism. A truncation of the DFT free energy functional up to three-point 
direct correlation functions has been used, and the dynamics derived from DDFT has been further simplified through 
lowest order approximations via a simple scale analysis to obtain the PFC equations, for both single-component 
and binary alloy systems. For the binary PFC model, the corresponding amplitude equations (both deterministic 
and stochastic) have been established via a hybrid multiple-scale approach, which describe large or "slow" scale 
dynamics of structural and compositional profiles based on the underlying crystalline state. Compared to other 
recent developments which have mainly focused on the evolution of complex structural amplitudes and concentration 
field, this work presents results that incorporate the new effects of mobility difference between alloy components, the 
coupling to zero-mode average atomic density, and also noise dynamics. Although the results of amplitude equations 
that we derive are for 2D hexagonal crystalline state, they can be extended to 3D bcc or fee structures by following 
a procedure similar to the one developed here and adopting the corresponding basic wavevectors (see also Ref . |3^] ) . 

This amplitude equation formalism for binary PFC has been applied to identifying the mechanisms and parameter 
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coupling during the process of surface segregation and alloy intermixing. Both liquid-solid and liquid-solid-solid 
epitaxial heterostructures have been examined, including morphological and compositional profiles. We find that the 
effect of concentration segregation on solid surface is controlled by material parameters such as the disparity of atomic 
size and mobility between different alloy components and misfit strain in solid layers. In the cases of nonzero solute 
expansion coefficient or unequal atomic mobility, an effect of composition overshooting around liquid-solid interface 
is obtained during strained layer growth, corresponding to vertical phase separation or segregation in the interface 
region. These results are consistent with recent experimental findings in heteroepitaxial systems, particularly the 
phenomenon of surface or interface segregation showing as the enrichment of one of the alloy species as compared to 
the bulk phase. This sample application of the amplitude equation formalism developed here has further illustrated 
the features and advantages of the PFC methodology, particularly in terms of modeling and understanding complex 
material phenomena involving spacial and temporal scales of experimental relevance. 
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Appendix A: Alternative derivations of binary PFC dynamics via DDFT 
1. Alternative derivation I 

In the following we provide an alternative derivation procedure for PFC dynamics, including two steps: i) directly 
use the original free energy functional ([T7]) and the DDFT equations to obtain the expressions of dp A{B)/ dt, and 
then ii) derive the dynamics of n and tp through Eq. p5|) . instead of using Eqs. (l2Tt and (|28)) as in Sec. IIIBl 

Define nA = {pA — pf')/pi and ub = [ps — pf)/Ph and hence the free energy functional ([T7|) can be rewritten as 
(using the expansion (|I^ 
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where Apf = pf / pi and Apf — pf / pi- From the DDFT equations (|I5p we can obtain the PFC equations for A & B 
components respectively, i.e.. 
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Note that n = ha + ns and ■0 



Kpf ~ pf) + PiiiT-A — nB)]/[pi{l + n)], and hence equivalent to Eq. ([25)1 we have 
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Substituting Eqs. (|X2|) and (Ml into Eq. (gH), and noting ua = (l + n)(l + '(/')/2- Ap^, ub = (l + n)(l- V')/2- Apf , 
A^A = Pi{Mi + M2), and Afs = pi{Mi — Af2), we can derive the binary PFC equations for n and t/j, which are exactly 
the same as Eqs. (gll), and ^U^. 



2. Alternative derivation II 



Another alternative derivation for PFC dynamics is to start with Eqs. as already derived in Sec. IIIBI 

Different from Sec. IIIBI in the formula ((28)) for Vi and I?2 if considering that the chemical potentials /i„ = SJ-/6n 
and fiN = 5T /5il) are slowly varying quantities and retaining terms up to the lowest order, we have 
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as used in the original PFC model [2^, [S^ ■ As in the previous work, the logarithm terms in the free energy functional 
(1211) are expanded in a power series, yielding (up to 4th order of n and ip) 
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where t = Bq ~ 1/2, w = 1 + /32, v = u = 1/3, and other parameters (such as B^, /3, K, and k) are defined in Eq. 

(1221). 

Substituting Eq. (jA6p into (jASp . we find (up to 3rd order) 
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As in Sec. IIIBI we can rescale the PFC equations (using the same length and time scales as well as the scales for n 
and ij} fields). Following the scale analysis discussed at the end of Sec. Ill Bl to lowest order approximation we can 
obtain the same simplified binary PFC equations given in Eqs. (|38|) and (p9|) . albeit with different forms of rescaled 
parameters 32 — 9ot/Bq and vi ~ 5o(l/2 + I3i)/Bq (other parameters g, wq, ao, and go are the same as those defined 
in Sec. Ill BI but with different value of w). 
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